R Notebook

library(forecast)
library(tseries)
library(lmtest)
library(Hmisc)
data <- read.csv("WAG_M.csv", sep=",", stringsAsFactors=F)
data <- head(data, -1)
names(data)[1] <- "Date"
names(data)[2] <- "Value"
data$Value <- as.numeric(data$Value)
data$Date <- as.Date(as.yearmon(data$Date, format="%Y-%m"))
tSeries <- ts(data = data$Value, start = as.numeric(c(format(data$Date[1], "%Y"), format(data$Date[1], "%m"))), freq = 12)
plot(tSeries, type="l", ylab="WAG_M", col="red")
grid()

trainSeries <- window(tSeries, end=c(2016,1))
testSeries  <- window(tSeries, start=c(2016,2))
D <- 24

STL-декомпозиция ряда:

plot(stl(tSeries, s.window="periodic"))

Оптимальное преобразование Бокса-Кокса и результат его применения:

par(mfrow=c(2,1))
plot(tSeries, ylab="Original series", xlab="", col="red")
grid()

LambdaOpt <- BoxCox.lambda(tSeries)
plot(BoxCox(tSeries, LambdaOpt), ylab="Transformed series", xlab="", col="red")
title(main=toString(round(LambdaOpt, 3)))

grid()

В данном случае преобразование имеет смысл использовать, так как оно хорошо стабилизирует дисперсию. Попробуем округлить параметр и взять \(\lambda=0.3\):

Результат практически такой же. Далее будем использовать \(\lambda=0.3\).

ARIMA

Автоматический подбор модели

Применим функцию auto.arima:

Series: tSeries 
ARIMA(2,1,2)(2,1,2)[12] 
Box Cox transformation: lambda= 0.2 

Coefficients:
         ar1      ar2      ma1     ma2     sar1    sar2    sma1     sma2
      1.7623  -0.8383  -1.6649  0.7626  -0.6909  0.0150  0.1017  -0.3917
s.e.  0.1902   0.1936   0.2308  0.2278   0.5095  0.1033  0.5038   0.2812

sigma^2 estimated as 0.006947:  log likelihood=300.88
AIC=-583.77   AICc=-583.1   BIC=-551.02

Предлагается модель ARIMA(0,1,0)(2,1,2)\(_{12}\). Посмотрим на её остатки:

Остаток симметричный относительно 0, похож на шум, значит хорошо отработало.

Гипотеза	Критерий	Результат проверки	Достигаемый уровень значимости
Нормальность	Шапиро-Уилка	отвергается	3.586421110^{-10} == 9.087472e-17
Несмещённость	Уилкоксона	отвергается	0.6895233 == 0.3847784
Стационарность	KPSS	отвергается	0.1 == 0.1

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

                      ME      RMSE       MAE        MPE     MAPE      MASE         ACF1 Theil's U
Training set -0.07669967  3.841494  2.553417 0.03242729 1.985713 0.1821145 -0.007336159        NA
Test set     13.01056720 14.883889 13.010567 5.48265396 5.482654 0.9279378  0.659315420 0.6510556

p-value greater than printed p-value

[1] 0.1

Ручной подбор модели

Исходный ряд нестационарен (p<0.01 == 0.01, критерий KPSS); сделаем сезонное дифференцирование:

Ряд стал стационарным (p<0.0105742==0.01057417, критерий KPSS).

Посмотрим на ACF и PACF полученного продифференцированного ряда:

На ACF значимы лаги 1 - 4, 44-47, на PACF — 1. Поищем с помощью auto.arima оптимальную модель полным перебором (stepwise=FALSE) с ограничениями d=1, D=1, max.p=5, max.q=3, max.P=0, max.Q=2, max.order=12:

Series: tSeries 
ARIMA(3,1,0)(0,1,1)[12] 
Box Cox transformation: lambda= 0.08800602 

Coefficients:
         ar1     ar2     ar3     sma1
      0.1611  0.0580  0.0949  -0.6465
s.e.  0.0598  0.0602  0.0605   0.0454

sigma^2 estimated as 0.00235:  log likelihood=453.83
AIC=-897.67   AICc=-897.45   BIC=-879.44

Была найдена более оптимальная по AICc модель — ARIMA(3,1,0)(0,1,1)\(_{12}\). Посмотрим на её остатки:

Отрежем начало ряда остатков и проанализируем их:

Достигаемые уровни значимости критерия Льюнга-Бокса для остатков:

Q-Q plot и гистограмма:

Гипотеза	Критерий	Результат проверки	Достигаемый уровень значимости
Нормальность	Шапиро-Уилка	отвергается	2.330884110^{-10} == 2.330884e-10
Несмещённость	Уилкоксона	не отвергается	0.5800223 == 0.5800223
Стационарность	KPSS	отвергается	0.1 == 0.1

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой: ARIMA(3,1,0)(0,1,1)\(_{12}\)

                      ME      RMSE       MAE          MPE     MAPE      MASE       ACF1 Theil's U
Training set -0.06702276  3.840327  2.541343 -0.007023565 1.971246 0.1812533 -0.1452584        NA
Test set     13.80180793 15.431423 13.801808  5.825917293 5.825917 0.9843706  0.6311357 0.6758536

Сравним остатки двух версий аримы, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:

    Diebold-Mariano Test

data:  resres.auto
DM = -0.11962, Forecast horizon = 1, Loss function power = 2, p-value = 0.9049
alternative hypothesis: two.sided

Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов.

В целом подобранная вручную модель проще, её остатки лучше, а ошибка на тесте меньше, так что остановимся на модели, подобранной вручную.